library("tidyverse")
library("firatheme")
library("prediction")
library("estimatr")
library("sandwich")
library("lmtest")
library("texreg")

# load data
load("~/df_replication.RData")

# Table SI 15 
model1 <- lm_robust(download ~ cue, data = df_replication, se_type="HC2")
summary(model1)
texreg(model1, digits = 3, include.ci = F)

# Table SI 16 
model2 <- lm_robust(download ~ issue, data = df_replication, se_type="HC2")
summary(model2)
texreg(model2, digits = 3, include.ci = F)

# Table SI 17 
model3 <- lm_robust(download ~ cue + issue + region + party_name, data = df_replication, se_type="HC2")
summary(model3)
texreg(model3, digits = 3, include.ci = F)

# Table SI 18 
model4 <- lm_robust(download ~ cue*issue + region + party_name, data = df_replication, se_type="HC2")
summary(model4)
texreg(model4, digits = 3, include.ci = F)


# Figure SI 15 - upper left

figure15.1 <- tibble(
  Treatment = c("Problem Indicators","Public Opinion", "Party Positions", "Media Reporting"),
  avg = c(round(coef(model1)[1]*100,1), round(coef(model1)[1]*100,1) + round(coef(model1)[2]*100,1),  round(coef(model1)[1]*100,1) + round(coef(model1)[3]*100,1),  round(coef(model1)[1]*100,1) + round(coef(model1)[4]*100,1)), 
  lower = c(round(coef(model1)[1]*100,1) - round(model1[["std.error"]][["(Intercept)"]]*100,4), round(coef(model1)[1]*100,1) + round(coef(model1)[2]*100,1) - round(model1[["std.error"]][["cuepublic_opinion"]]*100,4), round(coef(model1)[1]*100,1) + round(coef(model1)[3]*100,1) - round(model1[["std.error"]][["cuerival_parties"]]*100,4), round(coef(model1)[1]*100,1) + round(coef(model1)[4]*100,1) - round(model1[["std.error"]][["cuethe_media"]]*100,4)),
  upper = c(round(coef(model1)[1]*100,1) + round(model1[["std.error"]][["(Intercept)"]]*100,4), round(coef(model1)[1]*100,1) + round(coef(model1)[2]*100,1) + round(model1[["std.error"]][["cuepublic_opinion"]]*100,4), round(coef(model1)[1]*100,1) + round(coef(model1)[3]*100,1) + round(model1[["std.error"]][["cuerival_parties"]]*100,4), round(coef(model1)[1]*100,1) + round(coef(model1)[4]*100,1) + round(model1[["std.error"]][["cuethe_media"]]*100,4)))

ggplot(figure15.1, aes(x= reorder(Treatment,-avg),y=round(avg,1),fill=Treatment)) +
  geom_bar(stat="identity",position="dodge", color= "black", width=.25)+
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(.7))+
  geom_text(aes(y=upper, label=round(avg,3), size = 12), position=position_dodge(width=0.7), vjust=-.75)+
  labs(y="Download Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)")+
  scale_fill_manual(values=c("grey80","gray80", "grey80", "grey80")) + 
  scale_y_continuous(limits = c(0, 20)) + theme_fira() + 
  theme(text = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"))

# Figure SI 15 - upper right

figure15.2 <- tibble(
  Treatment = c("Climate Change","Public Schools"),
  avg = c(round(coef(model2)[1]*100,1), round(coef(model2)[1]*100,1) + round(coef(model2)[2]*100,1)), 
  lower = c(round(coef(model2)[1]*100,1) - round(model2[["std.error"]][["(Intercept)"]]*100,4), round(coef(model2)[1]*100,1) + round(coef(model2)[2]*100,1) - round(model2[["std.error"]][["issueschool"]]*100,4)),
  upper = c(round(coef(model2)[1]*100,1) + round(model2[["std.error"]][["(Intercept)"]]*100,4), round(coef(model2)[1]*100,1) + round(coef(model2)[2]*100,1) + round(model2[["std.error"]][["issueschool"]]*100,4)))

ggplot(figure15.2, aes(x= Treatment,y=avg, fill=Treatment)) +
  geom_bar(stat="identity",position="dodge", color= "black", width=.25)+
  geom_linerange(aes(ymin=lower, ymax=upper), position=position_dodge(.7))+
  geom_text(aes(y=upper, label=round(avg,3), size = 12), position=position_dodge(width=0.7), vjust=-.75)+
  labs(y="Download Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)", size = 4)+
  scale_fill_manual(values=c("grey80","gray50")) + 
  scale_y_continuous(limits = c(0, 20)) + theme_fira() + 
  theme(text = element_text(size=8), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"))


# Figure SI 15 - lower left 
estimates <- c(coef(model3)[2]*100,coef(model3)[3]*100,coef(model3)[4]*100,coef(model3)[5]*100)
lower95 <- c(confint(model3)[2,1]*100,confint(model3)[3,1]*100,confint(model3)[4,1]*100,confint(model3)[5,1]*100) 
upper95 <- c(confint(model3)[2,2]*100,confint(model3)[3,2]*100,confint(model3)[4,2]*100,confint(model3)[5,2]*100)  
lower90 <- c(confint(model3, level = 0.90)[2,1]*100,confint(model3, level = 0.90)[3,1]*100,confint(model3, level = 0.90)[4,1]*100,confint(model3, level = 0.90)[5,1]*100)
upper90 <- c(confint(model3, level = 0.90)[2,2]*100,confint(model3, level = 0.90)[3,2]*100,confint(model3, level = 0.90)[4,2]*100,confint(model3, level = 0.90)[5,2]*100) 
labels <- c("public opinion \n(vs. problem indicators)", "party position \n(vs. problem indicators)", "media reporting \n(vs. problem indicators)", "public schools \n(vs. climate change)")
type <- c("cue", "cue", "cue", "problem")
plotdf <- data.frame(estimates, lower95, upper95, lower90, upper90, labels, type)


ggplot(plotdf, aes(x=labels, y=estimates, group=labels, color = type)) + 
  geom_point(size=2, position=position_dodge(width=0.3)) +
  geom_linerange(aes(ymin=lower90, ymax=upper90), size=1.5, position=position_dodge(width=0.3)) + 
  geom_linerange(aes(ymin=lower95, ymax=upper95), size=0.5, position=position_dodge(width=0.3)) + 
  coord_flip() + theme_fira() + geom_hline(yintercept=0, linetype="dashed", 
                                           color = "gray50", size=0.5) +
  scale_color_manual(values=c("gray50", "gray50")) +
  theme(legend.position= "none") + labs(x="", y="Δ Percentage points", title="", caption = "Note: Error bars represent 90% and 95% Confidence Intervals") +
  ylim(-10, 10) + facet_wrap(~type, ncol = 1, scales = "free_y") +
  theme(axis.title = element_text(size=6), legend.position= "none", plot.caption = element_text(size = 6, color = "gray40", face = "italic"))

# Figure SI 15 - lower right

a <- prediction(model3, at = list(issue = "school", cue = "problem_indicators"))
b <- prediction(model3, at = list(issue = "climate", cue = "problem_indicators"))
c <- prediction(model3, at = list(issue = "school", cue = "public_opinion"))
d <- prediction(model3, at = list(issue = "climate", cue = "public_opinion"))
e <- prediction(model3, at = list(issue = "school", cue = "rival_parties"))
f <- prediction(model3, at = list(issue = "climate", cue = "rival_parties"))
g <- prediction(model3, at = list(issue = "school", cue = "the_media"))
h <- prediction(model3, at = list(issue = "climate", cue = "the_media"))


figure15.4 <- tibble(
  Treatment = c("Problem Indicators","Problem Indicators","Public Opinion", "Public Opinion", "Party Positions", "Party Positions", "Media Reporting", "Media Reporting"),
  Problem = c("Public Schools","Climate Change","Public Schools", "Climate Change", "Public Schools", "Climate Change", "Public Schools", "Climate Change"),
  avg = c(mean(a$fitted)*100, mean(b$fitted)*100, mean(c$fitted)*100, mean(d$fitted)*100, mean(e$fitted)*100, mean(f$fitted)*100, mean(g$fitted)*100, mean(h$fitted)*100), 
  lower = c(mean(a$fitted)*100 - mean(a$se.fitted)*100, mean(b$fitted)*100 - mean(b$se.fitted)*100, mean(c$fitted)*100 - mean(c$se.fitted)*100, mean(d$fitted)*100 - mean(d$se.fitted)*100, mean(e$fitted)*100 - mean(e$se.fitted)*100, mean(f$fitted)*100 - mean(f$se.fitted)*100, mean(g$fitted)*100 - mean(g$se.fitted)*100, mean(h$fitted)*100 - mean(h$se.fitted)*100),
  upper = c(mean(a$fitted)*100 + mean(a$se.fitted)*100, mean(b$fitted)*100 + mean(b$se.fitted)*100, mean(c$fitted)*100 + mean(c$se.fitted)*100, mean(d$fitted)*100 + mean(d$se.fitted)*100, mean(e$fitted)*100 + mean(e$se.fitted)*100, mean(f$fitted)*100 + mean(f$se.fitted)*100, mean(g$fitted)*100 + mean(g$se.fitted)*100, mean(h$fitted)*100 + mean(h$se.fitted)*100))

ggplot(figure15.4, aes(x= reorder(Treatment,-avg), y=avg, fill=Problem)) +
  geom_bar(stat="identity", position="dodge", color= "black", width=.25) +
  geom_linerange(aes(ymin=lower, ymax=upper), position= position_dodge(.25)) +
  geom_text(aes(y=upper, label=round(avg,1)), size = 3, position=position_dodge(width=0.7), vjust=-1.75)+
  labs(y="Download Rate (%)", x="", title="", caption = "Note: Error bars represent ± 1 standard error (HC2)", size = 4)+
  scale_fill_manual(values=c("grey80","gray50", "grey80","gray50", "grey30","gray50", "grey80","gray50")) + 
  scale_y_continuous(limits = c(0, 20)) + theme_fira() + 
  theme(text = element_text(size=8), legend.position= "bottom", plot.caption = element_text(size = 6, color = "gray40", face = "italic"))
